pSERG CIRCADIAN DISTRIBUTION OF STATUS EPILEPTICUS

---

This analysis evaluates the circadian distribution of refractory status epilepticus episodes

Install packages needed for analysis

# install.packages("plyr")
library(plyr)

# install.packages("gdata")
library(gdata)

## gdata: Unable to locate valid perl interpreter
## gdata: 
## gdata: read.xls() will be unable to read Excel XLS and XLSX files
## gdata: unless the 'perl=' argument is used to specify the location
## gdata: of a valid perl intrpreter.
## gdata: 
## gdata: (To avoid display of this message in the future, please
## gdata: ensure perl is installed and available on the executable
## gdata: search path.)

## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLX' (Excel 97-2004) files.

## 

## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLSX' (Excel 2007+) files.

## 

## gdata: Run the function 'installXLSXsupport()'
## gdata: to automatically download and install the perl
## gdata: libaries needed to support Excel XLS and XLSX formats.

## 
## Attaching package: 'gdata'

## The following object is masked from 'package:stats':
## 
##     nobs

## The following object is masked from 'package:utils':
## 
##     object.size

## The following object is masked from 'package:base':
## 
##     startsWith

# install.packages("gmodels")
library(gmodels)

# install.packages("cosinor")
library(cosinor)

## Warning: package 'cosinor' was built under R version 3.5.3

# install.packages("ggplot2")
library(ggplot2)

# install.packages("plotly")
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following objects are masked from 'package:plyr':
## 
##     arrange, mutate, rename, summarise

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

Load the database

## Load the data
pSERGcircadian <- read.csv("F:\\PROJECTS\\pSERG Circadian\\Data\\pSERG.csv")

Data cleaning and preparation of variables for analysis

## Keep patients only if they are a refractory SE case (not a control)
pSERGcircadian <- pSERGcircadian[which(pSERGcircadian$SE_GROUP == 'refractory_case'), ]
pSERGcircadian$SE_GROUP <- droplevels(pSERGcircadian$SE_GROUP)

## Transform date of status epilepticus into date format
pSERGcircadian$DATESEIZURE <- as.Date(pSERGcircadian$DATESEIZURE, format = "%m/%d/%Y")

## Order by date of status epilepticus
pSERGcircadian <- pSERGcircadian[order(pSERGcircadian$PATIENT_LABEL, pSERGcircadian$DATESEIZURE), ]

## Delete duplicate episodes from the same patient
pSERGcircadian <- pSERGcircadian[!duplicated(pSERGcircadian$PATIENT_LABEL), ]


#################Create variable hour of SE######################
## Transform time of hour to numeric and delete cases with unknown hour of SE
pSERGcircadian$TIMESEIZURE_HOURS <- as.numeric(as.character(pSERGcircadian$TIMESEIZURE_HOURS))

## Warning: NAs introduced by coercion

pSERGcircadian$TIMESEIZURE_HOURS <- ifelse(pSERGcircadian$TIMESEIZURE_HOURS>99, pSERGcircadian$TIMESEIZURE_HOURS%/%100, pSERGcircadian$TIMESEIZURE_HOURS)

## Delete patients with no numeric values for hour of seizure onset or no information
pSERGcircadian <- pSERGcircadian[!is.na(pSERGcircadian$TIMESEIZURE_HOURS), ]

## Rename variable
pSERGcircadian$hourofSE <- pSERGcircadian$TIMESEIZURE_HOURS

## Rename 0 to 24
pSERGcircadian$hourofSE[pSERGcircadian$hourofSE == 0] <- 24

## Exact hour
pSERGcircadian$TIMESEIZURE_MINUTES <- as.numeric(as.character(pSERGcircadian$TIMESEIZURE_MINUTES))

## Warning: NAs introduced by coercion

pSERGcircadian$time <- pSERGcircadian$TIMESEIZURE_HOURS + (pSERGcircadian$TIMESEIZURE_MINUTES / 60)


#######################Transform age into numeric and delete patients with unknown age####################
# Transform age into numeric
pSERGcircadian$G_T_STTS_PLPTCS_EPISODE_MONTHS <- as.numeric(as.character(pSERGcircadian$G_T_STTS_PLPTCS_EPISODE_MONTHS))

## Warning: NAs introduced by coercion

pSERGcircadian$G_T_STTS_PLPTCUS_EPISODE_YEARS <- as.numeric(as.character(pSERGcircadian$G_T_STTS_PLPTCUS_EPISODE_YEARS))

## Warning: NAs introduced by coercion

pSERGcircadian$ageyears <- pSERGcircadian$G_T_STTS_PLPTCUS_EPISODE_YEARS + (pSERGcircadian$G_T_STTS_PLPTCS_EPISODE_MONTHS/12)
pSERGcircadian <- pSERGcircadian[complete.cases(pSERGcircadian[ , "ageyears"]), ]


##########################Delete patients with unknown sex################################################
pSERGcircadian <- pSERGcircadian[which(pSERGcircadian$SEX == "male" | pSERGcircadian$SEX == "female"), ]


#######################Transform variables for analysis###############
## Create variable delay
pSERGcircadian$delay[grepl("delay", pSERGcircadian$PAST)] <- 1
pSERGcircadian$delay[!grepl("delay", pSERGcircadian$PAST)] <- 0

## Create variable palsy
pSERGcircadian$palsy[grepl("palsy", pSERGcircadian$PAST)] <- 1
pSERGcircadian$palsy[!grepl("palsy", pSERGcircadian$PAST)] <- 0
 
## Create variable febrile
pSERGcircadian$febrile[grepl("febrile", pSERGcircadian$PAST)] <- 1
pSERGcircadian$febrile[!grepl("febrile", pSERGcircadian$PAST)] <- 0

## Create variable prior epilepsy
pSERGcircadian$priorepilepsy[grepl("epi", pSERGcircadian$PAST)] <- 1
pSERGcircadian$priorepilepsy[!grepl("epi", pSERGcircadian$PAST)] <- 0

## Create variable status
pSERGcircadian$status[grepl("status", pSERGcircadian$PAST)] <- 1
pSERGcircadian$status[!grepl("status", pSERGcircadian$PAST)] <- 0

## Create variable none
pSERGcircadian$none[grepl("none", pSERGcircadian$PAST)] <- 1
pSERGcircadian$none[!grepl("none", pSERGcircadian$PAST)] <- 0

Demographics

############################DEMOGRAPHICS#######################
## Gender
CrossTable(pSERGcircadian$SEX)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##           |    female |      male | 
##           |-----------|-----------|
##           |       155 |       213 | 
##           |     0.421 |     0.579 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Age
nobs(pSERGcircadian$ageyears)

## [1] 368

summary(pSERGcircadian$ageyears)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.08333  1.32146  4.16667  5.94850  9.68540 20.74167

sd(pSERGcircadian$ageyears, na.rm = TRUE)

## [1] 5.27164

# Patients 18 year-old or older
pSERGcircadian$eighteenormore <- as.factor(ifelse(pSERGcircadian$ageyears >= 18, 1, 0))
CrossTable(pSERGcircadian$eighteenormore)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       361 |         7 | 
##           |     0.981 |     0.019 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Race
CrossTable(pSERGcircadian$RACE)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##                                     |       american_indian_alaska_native |                              arabic |                               asian |           black_or_african_american | native_hawaiian_or_pacific_islander | 
##                                     |-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|
##                                     |                                   1 |                                  10 |                                  11 |                                  71 |                                   1 | 
##                                     |                               0.003 |                               0.027 |                               0.030 |                               0.193 |                               0.003 | 
##                                     |-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|
## 
## 
##                                     |                        not_reported |                             unknown |                               white | 
##                                     |-------------------------------------|-------------------------------------|-------------------------------------|
##                                     |                                  15 |                                  22 |                                 237 | 
##                                     |                               0.041 |                               0.060 |                               0.644 | 
##                                     |-------------------------------------|-------------------------------------|-------------------------------------|
## 
## 
## 
## 

## Ethnicity
CrossTable(pSERGcircadian$ETHNICITY)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##                        |     hispanic_or_latino | not_hispanic_or_latino |           not_reported |                unknown | 
##                        |------------------------|------------------------|------------------------|------------------------|
##                        |                     66 |                    269 |                     21 |                     12 | 
##                        |                  0.179 |                  0.731 |                  0.057 |                  0.033 | 
##                        |------------------------|------------------------|------------------------|------------------------|
## 
## 
## 
## 

## Delay
CrossTable(pSERGcircadian$delay)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       182 |       186 | 
##           |     0.495 |     0.505 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Palsy
CrossTable(pSERGcircadian$palsy)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       329 |        39 | 
##           |     0.894 |     0.106 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Febrile
CrossTable(pSERGcircadian$febrile)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       328 |        40 | 
##           |     0.891 |     0.109 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Prior epilepsy
CrossTable(pSERGcircadian$priorepilepsy)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       186 |       182 | 
##           |     0.505 |     0.495 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Status
CrossTable(pSERGcircadian$status)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       299 |        69 | 
##           |     0.812 |     0.188 | 
##           |-----------|-----------|
## 
## 
## 
## 

## None
CrossTable(pSERGcircadian$none)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       248 |       120 | 
##           |     0.674 |     0.326 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Transform CONVULSIVEDURATION into numeric
pSERGcircadian$CONVULSIVEDURATION <- as.numeric(as.character(pSERGcircadian$CONVULSIVEDURATION))

## Warning: NAs introduced by coercion

pSERGcircadian$CONVULSIVEmin <- pSERGcircadian$CONVULSIVEDURATION
pSERGcircadian$CONVULSIVEhr <- pSERGcircadian$CONVULSIVEDURATION*60
pSERGcircadian$convulsivedurationmin <- ifelse(pSERGcircadian$CONVULSIVEDURATIONUNITS == "min", pSERGcircadian$CONVULSIVEmin, pSERGcircadian$CONVULSIVEhr)

## Duration of SE
nobs(pSERGcircadian$convulsivedurationmin)

## [1] 355

summary(pSERGcircadian$convulsivedurationmin)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##       0      60     122    1925     283  172800      13

sd(pSERGcircadian$convulsivedurationmin, na.rm = TRUE)

## [1] 12308.02

## Hospital onset
CrossTable(pSERGcircadian$HOSPITALONSET)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##           |        no |       yes | 
##           |-----------|-----------|
##           |       258 |       110 | 
##           |     0.701 |     0.299 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Type of status epilepticus
CrossTable(pSERGcircadian$TYPESTATUS)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  368 
## 
##  
##              |   continuous | intermittent | 
##              |--------------|--------------|
##              |          116 |          252 | 
##              |        0.315 |        0.685 | 
##              |--------------|--------------|
## 
## 
## 
## 

Figures of the distribution of status epilepticus episodes

## Create dummy variable for events of SE
pSERGcircadian$events <- 1


################Every hour###############
## Draw the number of SE episodes in a plot with stacked units
distributionhour <- ggplot(pSERGcircadian, aes(x = hourofSE, y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") + 
  scale_x_discrete(limits = c(1,24)) +
  theme(panel.background = element_rect(fill = "white"), 
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15, 20, 25)) +
  labs(x= "Time of day (hour)", y= "Number of rSE episodes")

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

distributionhour

# Interactive graph
ggplotly(distributionhour)

################Every 2 hours###############
# Create variable every 2 hours
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 1] <- "1-2"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 2] <- "1-2"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 3] <- "3-4"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 4] <- "3-4"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 5] <- "5-6"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 6] <- "5-6"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 7] <- "7-8"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 8] <- "7-8"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 9] <- "9-10"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 10] <- "9-10"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 11] <- "11-12"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 12] <- "11-12"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 13] <- "13-14"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 14] <- "13-14"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 15] <- "15-16"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 16] <- "15-16"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 17] <- "17-18"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 18] <- "17-18"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 19] <- "19-20"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 20] <- "19-20"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 21] <- "21-22"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 22] <- "21-22"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 23] <- "23-24"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 24] <- "23-24"

# Give this variable an order
pSERGcircadian$groups2 <- factor(pSERGcircadian$groups2, c("1-2", "3-4","5-6", "7-8","9-10", "11-12","13-14", "15-16","17-18", "19-20","21-22", "23-24"))

## Draw the number of SE episodes in a plot with stacked units
distribution2hours <- ggplot(pSERGcircadian[!is.na(pSERGcircadian$groups2), ], aes(x = groups2, y = events))+geom_bar(stat ='identity', color = "#0076c0", fill = "#67771a") + 
  theme(panel.background = element_rect(fill = "white"), 
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_discrete() + scale_y_continuous(breaks = c(10, 20, 30, 40)) +
  labs(x= "Time of day (hour)", y= "Number of rSE episodes")
distribution2hours

# Interactive graph
ggplotly(distribution2hours)

################Every 3 hours###############
# Create variable every 3 hours
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 1] <- "1-3"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 2] <- "1-3"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 3] <- "1-3"

pSERGcircadian$groups3[pSERGcircadian$hourofSE == 4] <- "4-6"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 5] <- "4-6"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 6] <- "4-6"

pSERGcircadian$groups3[pSERGcircadian$hourofSE == 7] <- "7-9"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 8] <- "7-9"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 9] <- "7-9"

pSERGcircadian$groups3[pSERGcircadian$hourofSE == 10] <- "10-12"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 11] <- "10-12"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 12] <- "10-12"

pSERGcircadian$groups3[pSERGcircadian$hourofSE == 13] <- "13-15"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 14] <- "13-15"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 15] <- "13-15"

pSERGcircadian$groups3[pSERGcircadian$hourofSE== 16] <- "16-18"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 17] <- "16-18"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 18] <- "16-18"

pSERGcircadian$groups3[pSERGcircadian$hourofSE== 19] <- "19-21"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 20] <- "19-21"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 21] <- "19-21"
 
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 22] <- "22-24"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 23] <- "22-24"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 24] <- "22-24"

# Give this variable an order
pSERGcircadian$groups3 <- factor(pSERGcircadian$groups3, c("1-3", "4-6", "7-9", "10-12", "13-15", "16-18", "19-21", "22-24"))

## Draw the number of SE in a plot
distribution3hours <- ggplot(pSERGcircadian[!is.na(pSERGcircadian$groups3), ], aes(x = groups3, y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") + 
  theme(panel.background=element_rect(fill = "white"), 
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_discrete() + scale_y_continuous(breaks = c(10, 20, 30, 40, 50)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes")
distribution3hours

# Interactive graph
ggplotly(distribution3hours)

################Every 4 hours###############
#Create variable every 4 hours
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 1] <- "1-4"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 2] <- "1-4"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 3] <- "1-4"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 4] <- "1-4"

pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 5] <- "5-8"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 6] <- "5-8"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 7] <- "5-8"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 8] <- "5-8"

pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 9] <- "9-12"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 10] <- "9-12"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 11] <- "9-12"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 12] <- "9-12"

pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 13] <- "13-16"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 14] <- "13-16"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 15] <- "13-16"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 16] <- "13-16"

pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 17] <- "17-20"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 18] <- "17-20"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 19] <- "17-20"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 20] <- "17-20"

pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 21] <- "21-24"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 22] <- "21-24"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 23] <- "21-24"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 24] <- "21-24"

#Give this variable an order
pSERGcircadian$groupsfour <- factor(pSERGcircadian$groupsfour, c("1-4","5-8","9-12","13-16","17-20","21-24"))

##Draw the number of SE in a plot
distribution4hours <- ggplot(pSERGcircadian[!is.na(pSERGcircadian$groupsfour), ], aes(x = groupsfour, y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") + 
  theme(panel.background = element_rect(fill = "white"), 
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_discrete() + scale_y_continuous(breaks = c(10, 20, 30, 40, 50, 60, 70)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes")
distribution4hours

# Interactive graph
ggplotly(distribution4hours)

Primary outcome

# Test (Kolmogorov-Smirnov) if the distribution of SE during the day is consistent with a uniform distribution
ks.test(pSERGcircadian$hourofSE, "punif", 0, 24)

## Warning in ks.test(pSERGcircadian$hourofSE, "punif", 0, 24): ties should
## not be present for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  pSERGcircadian$hourofSE
## D = 0.088768, p-value = 0.006058
## alternative hypothesis: two-sided

####################Distribution of SE episodes through the day###########################
# Distribution of SE episodes if they were evenly distributed
nobs(pSERGcircadian$events)

## [1] 368

nobs(pSERGcircadian$events)/24

## [1] 15.33333

#####################################################################################
# Create another database with the circadian data
occurrences <- count(pSERGcircadian, "hourofSE")
# Make it a dataframe
circadiandata <- data.frame(occurrences)
circadiandata$numberSE <- circadiandata$freq
circadiandata <- circadiandata[ , c("hourofSE", "numberSE")]
#####################################################################################



# Cosinor analysis
SEfit <- cosinor.lm(numberSE ~ time(hourofSE), data = circadiandata, period = 24)
summary(SEfit)

## Raw model coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  15.3333         0.7472  13.8688  16.7978  0.0000
## rrr          -3.1912         1.0567  -5.2623  -1.1200  0.0025
## sss           0.3806         1.0567  -1.6905   2.4517  0.7187
## 
## ***********************
## 
## Transformed coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  15.3333         0.7472  13.8688  16.7978  0.0000
## amp           3.2138         1.0567   1.1427   5.2849  0.0024
## acr          -0.1187         0.3288  -0.7631   0.5257  0.7181

# Draw the fitting curve 1h
## Draw the number of SE in a plot
circadianfigure <- ggplot.cosinor.lm(SEfit) + geom_bar(data = pSERGcircadian, aes(x = hourofSE, y = events), stat = 'identity', color = "#0076c0", fill = "#67771a", alpha = 0.8) + 
  scale_x_discrete(limits = c(1,24)) +
  theme(panel.background = element_rect(fill = "white"), 
        axis.text=element_text(size = 18, color = "black", face = "bold"),
        axis.title=element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15, 20, 25)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
  geom_segment(aes(x = 0, xend = 24, y = SEfit$coefficients[1], yend = SEfit$coefficients[1])) +
  geom_segment(aes(x = 11.45, xend = 11.45, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] + SEfit$coefficients[2])), size=1.2) +
  geom_segment(aes(x = 23.45, xend = 23.45, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] - SEfit$coefficients[2])), size=1.2)

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

circadianfigure

# Interactive graph
ggplotly(circadianfigure)

#Draw the data and the fitting curve
circadianonly <- ggplot.cosinor.lm(SEfit) +
  geom_point(data = circadiandata, aes(hourofSE, numberSE), color = "#67771a", size = 4) +
  scale_x_discrete(limits = c(1,24)) + 
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
        axis.line = element_line(color= "black", size = 1.2)) +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
  geom_segment(aes(x = 0, xend = 24, y = SEfit$coefficients[1], yend = SEfit$coefficients[1])) +
  geom_segment(aes(x = 10.55, xend = 10.55, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] + SEfit$coefficients[2]))) +
  geom_segment(aes(x = 22.55, xend = 22.55, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] - SEfit$coefficients[2])))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

circadianonly

# Interactive graph
ggplotly(circadianonly)

Comparison patients with and without a neurologic history

## Proportion of patients with some neurological history per time bin
neurologichistory1h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$none)[2])
neurologichistory2h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$none)[2])
neurologichistory3h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$none)[2])
neurologichistory4h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$none)[2])
neurologichistory5h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$none)[2])
neurologichistory6h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$none)[2])
neurologichistory7h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$none)[2])
neurologichistory8h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$none)[2])
neurologichistory9h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$none)[2])
neurologichistory10h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$none)[2])
neurologichistory11h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$none)[2])
neurologichistory12h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$none)[2])
neurologichistory13h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$none)[2])
neurologichistory14h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$none)[2])
neurologichistory15h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$none)[2])
neurologichistory16h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$none)[2])
neurologichistory17h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$none)[2])
neurologichistory18h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$none)[2])
neurologichistory19h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$none)[2])
neurologichistory20h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$none)[2])
neurologichistory21h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$none)[2])
neurologichistory22h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$none)[2])
neurologichistory23h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$none)[2])
neurologichistory24h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$none)[2])

circadiandata$neurologic <- c(neurologichistory1h, neurologichistory2h, neurologichistory3h, neurologichistory4h, neurologichistory5h, neurologichistory6h, neurologichistory7h, neurologichistory8h, neurologichistory9h, neurologichistory10h, neurologichistory11h, neurologichistory12h, neurologichistory13h, neurologichistory14h, neurologichistory15h, neurologichistory16h, neurologichistory17h, neurologichistory18h, neurologichistory19h, neurologichistory20h, neurologichistory21h, neurologichistory22h, neurologichistory23h, neurologichistory24h)
circadiandata$neurologic <- ifelse(is.na(circadiandata$neurologic) == TRUE, 0, circadiandata$neurologic)
circadiandata$numberneurologic <- circadiandata$numberSE * circadiandata$neurologic


## Cosinor analysis
SEfitcomparisonneurologic <- cosinor.lm(numberSE ~ time(hourofSE) + numberneurologic + amp.acro(numberneurologic), data = circadiandata, period = 24)
summary(SEfitcomparisonneurologic)

## Raw model coefficients:
##                      estimate standard.error lower.CI upper.CI p.value
## (Intercept)            5.3261         1.6303   2.1307   8.5215  0.0011
## numberneurologic       0.9541         0.1597   0.6412   1.2671  0.0000
## rrr                    0.0079         1.7118  -3.3471   3.3630  0.9963
## sss                   -1.4530         2.9021  -7.1411   4.2351  0.6166
## numberneurologic:rrr  -0.2224         0.1635  -0.5427   0.0980  0.1737
## numberneurologic:sss   0.0512         0.2722  -0.4823   0.5846  0.8509
## 
## ***********************
## 
## Transformed coefficients:
##                            estimate standard.error lower.CI upper.CI
## (Intercept)                  5.3261         1.6303   2.1307   8.5215
## [numberneurologic = 1]       0.9541         0.1597   0.6412   1.2671
## amp                          1.4530         2.9012  -4.2332   7.1392
## [numberneurologic = 1]:amp   1.4181         2.6385  -3.7533   6.5896
## acr                         -1.5653         1.1793  -3.8766   0.7460
## [numberneurologic = 1]:acr   1.4190         1.0970  -0.7311   3.5691
##                            p.value
## (Intercept)                 0.0011
## [numberneurologic = 1]      0.0000
## amp                         0.6165
## [numberneurologic = 1]:amp  0.5909
## acr                         0.1844
## [numberneurologic = 1]:acr  0.1958

test_cosinor(SEfitcomparisonneurologic, "numberneurologic", param = "amp")

## Global test: 
## Statistic: 
## [1] 0.01
## 
## 
##  P-value: 
## [1] 0.9234
## 
##  Individual tests: 
## Statistic: 
## [1] -0.1
## 
## 
##  P-value: 
## [1] 0.9234
## 
##  Estimate and confidence interval[1] "-0.03 (-0.75 to 0.68)"

test_cosinor(SEfitcomparisonneurologic, "numberneurologic", param = "acr")

## Global test: 
## Statistic: 
## [1] 86.37
## 
## 
##  P-value: 
## [1] 0
## 
##  Individual tests: 
## Statistic: 
## [1] 9.29
## 
## 
##  P-value: 
## [1] 0
## 
##  Estimate and confidence interval[1] "2.98 (2.35 to 3.61)"

##Graph
neurologichistoryfigure <- ggplot.cosinor.lm(SEfitcomparisonneurologic, x_str = "numberneurologic") +
  geom_bar(data = pSERGcircadian[pSERGcircadian$none == 0, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#5698a3", fill = "#5698a3", alpha = 0.2) + 
  geom_bar(data = pSERGcircadian[pSERGcircadian$none == 1, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#a30234", fill = "#a30234", alpha = 0.2) +
    scale_x_discrete(limits = c(1,24)) +
  theme(panel.background=element_rect(fill = "white"), 
        axis.text=element_text(size = 18, color = "black", face = "bold"),
        axis.title=element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length=unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
        legend.position = "none") +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15, 20)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonneurologic$coefficients[1], yend = SEfitcomparisonneurologic$coefficients[1]), color = "#a30234") +
  geom_segment(aes(x = 18.2, xend = 18.2, y = SEfitcomparisonneurologic$coefficients[1], yend = (SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[3])), color = "#a30234", size=1.2) +
  geom_segment(aes(x = 6.2, xend = 6.2, y = SEfitcomparisonneurologic$coefficients[1], yend = (SEfitcomparisonneurologic$coefficients[1] - SEfitcomparisonneurologic$coefficients[3])), color = "#a30234", size=1.2) +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2], yend = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2]), color = "#5698a3") +
  geom_segment(aes(x = 17.4, xend = 17.4, y = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2], yend = (SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2] + SEfitcomparisonneurologic$coefficients[4])), color = "#5698a3", size=1.2) + geom_segment(aes(x = 5.4, xend = 5.4, y = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2], yend = (SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2] - SEfitcomparisonneurologic$coefficients[4])), color = "#5698a3", size=1.2)

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neurologichistoryfigure

# Interactive graph
ggplotly(neurologichistoryfigure)

Comparison patients with and without a history of epilepsy

## Proportion of patients with a history of epilepsy per time bin
epilepsyhistory1h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$priorepilepsy)[2])
epilepsyhistory2h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$priorepilepsy)[2])
epilepsyhistory3h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$priorepilepsy)[2])
epilepsyhistory4h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$priorepilepsy)[2])
epilepsyhistory5h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$priorepilepsy)[2])
epilepsyhistory6h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$priorepilepsy)[2])
epilepsyhistory7h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$priorepilepsy)[2])
epilepsyhistory8h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$priorepilepsy)[2])
epilepsyhistory9h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$priorepilepsy)[2])
epilepsyhistory10h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$priorepilepsy)[2])
epilepsyhistory11h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$priorepilepsy)[2])
epilepsyhistory12h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$priorepilepsy)[2])
epilepsyhistory13h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$priorepilepsy)[2])
epilepsyhistory14h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$priorepilepsy)[2])
epilepsyhistory15h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$priorepilepsy)[2])
epilepsyhistory16h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$priorepilepsy)[2])
epilepsyhistory17h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$priorepilepsy)[2])
epilepsyhistory18h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$priorepilepsy)[2])
epilepsyhistory19h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$priorepilepsy)[2])
epilepsyhistory20h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$priorepilepsy)[2])
epilepsyhistory21h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$priorepilepsy)[2])
epilepsyhistory22h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$priorepilepsy)[2])
epilepsyhistory23h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$priorepilepsy)[2])
epilepsyhistory24h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$priorepilepsy)[2])

circadiandata$epilepsy <- c(epilepsyhistory1h, epilepsyhistory2h, epilepsyhistory3h, epilepsyhistory4h, epilepsyhistory5h, epilepsyhistory6h, epilepsyhistory7h, epilepsyhistory8h, epilepsyhistory9h, epilepsyhistory10h, epilepsyhistory11h, epilepsyhistory12h, epilepsyhistory13h, epilepsyhistory14h, epilepsyhistory15h, epilepsyhistory16h, epilepsyhistory17h, epilepsyhistory18h, epilepsyhistory19h, epilepsyhistory20h, epilepsyhistory21h, epilepsyhistory22h, epilepsyhistory23h, epilepsyhistory24h)
circadiandata$epilepsy <- ifelse(is.na(circadiandata$epilepsy) == TRUE, 0, circadiandata$epilepsy)
circadiandata$numberepilepsy <- circadiandata$numberSE * circadiandata$epilepsy

## Cosinor analysis
SEfitcomparisonepilepsy <- cosinor.lm(numberSE ~ time(hourofSE) + numberepilepsy + amp.acro(numberepilepsy), data = circadiandata, period = 24)
summary(SEfitcomparisonepilepsy)

## Raw model coefficients:
##                    estimate standard.error lower.CI upper.CI p.value
## (Intercept)          7.4248         1.4394   4.6037  10.2460  0.0000
## numberepilepsy       1.0437         0.2000   0.6518   1.4357  0.0000
## rrr                 -1.3100         1.4850  -4.2206   1.6006  0.3777
## sss                 -0.6892         2.7395  -6.0585   4.6800  0.8014
## numberepilepsy:rrr  -0.1829         0.1869  -0.5492   0.1834  0.3277
## numberepilepsy:sss  -0.0757         0.3557  -0.7728   0.6213  0.8314
## 
## ***********************
## 
## Transformed coefficients:
##                          estimate standard.error lower.CI upper.CI p.value
## (Intercept)                7.4248         1.4394   4.6037  10.2460  0.0000
## [numberepilepsy = 1]       1.0437         0.2000   0.6518   1.4357  0.0000
## amp                        1.4803         1.7568  -1.9630   4.9235  0.3994
## [numberepilepsy = 1]:amp   1.6775         1.5324  -1.3260   4.6810  0.2737
## acr                        0.4843         1.7386  -2.9233   3.8920  0.7806
## [numberepilepsy = 1]:acr   0.4735         1.3506  -2.1735   3.1206  0.7259

test_cosinor(SEfitcomparisonepilepsy, "numberepilepsy", param = "amp")

## Global test: 
## Statistic: 
## [1] 0.66
## 
## 
##  P-value: 
## [1] 0.4182
## 
##  Individual tests: 
## Statistic: 
## [1] 0.81
## 
## 
##  P-value: 
## [1] 0.4182
## 
##  Estimate and confidence interval[1] "0.2 (-0.28 to 0.67)"

test_cosinor(SEfitcomparisonepilepsy, "numberepilepsy", param = "acr")

## Global test: 
## Statistic: 
## [1] 0
## 
## 
##  P-value: 
## [1] 0.978
## 
##  Individual tests: 
## Statistic: 
## [1] -0.03
## 
## 
##  P-value: 
## [1] 0.978
## 
##  Estimate and confidence interval[1] "-0.01 (-0.78 to 0.76)"

##Graph
epilepsyhistoryfigure <- ggplot.cosinor.lm(SEfitcomparisonepilepsy, x_str = "numberepilepsy") +
  geom_bar(data = pSERGcircadian[pSERGcircadian$priorepilepsy == 1, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#5698a3", fill = "#5698a3", alpha = 0.2) + 
  geom_bar(data = pSERGcircadian[pSERGcircadian$priorepilepsy == 0, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#a30234", fill = "#a30234", alpha = 0.2) +
    scale_x_discrete(limits = c(1,24)) +
  theme(panel.background=element_rect(fill = "white"), 
        axis.text=element_text(size = 18, color = "black", face = "bold"),
        axis.title=element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length=unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
        legend.position = "none") +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonepilepsy$coefficients[1], yend = SEfitcomparisonepilepsy$coefficients[1]), color = "#a30234") +
  geom_segment(aes(x = 13.9, xend = 13.9, y = SEfitcomparisonepilepsy$coefficients[1], yend = (SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[3])), color = "#a30234", size=1.2) +
  geom_segment(aes(x = 1.9, xend = 1.9, y = SEfitcomparisonepilepsy$coefficients[1], yend = (SEfitcomparisonepilepsy$coefficients[1] - SEfitcomparisonepilepsy$coefficients[3])), color = "#a30234", size=1.2) +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2], yend = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2]), color = "#5698a3") +
  geom_segment(aes(x = 13.8, xend = 13.8, y = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2], yend = (SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2] + SEfitcomparisonepilepsy$coefficients[4])), color = "#5698a3", size=1.2) + geom_segment(aes(x = 1.8, xend = 1.8, y = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2], yend = (SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2] - SEfitcomparisonepilepsy$coefficients[4])), color = "#5698a3", size=1.2)

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

epilepsyhistoryfigure

# Interactive graph
ggplotly(epilepsyhistoryfigure)

Comparison patients with SE onset in and out of the hospital

## Proportion of patients with onset of SE out of the hospital
outofhospitalonset1h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$HOSPITALONSET)[2])
outofhospitalonset2h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$HOSPITALONSET)[2])
outofhospitalonset3h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$HOSPITALONSET)[2])
outofhospitalonset4h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$HOSPITALONSET)[2])
outofhospitalonset5h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$HOSPITALONSET)[2])
outofhospitalonset6h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$HOSPITALONSET)[2])
outofhospitalonset7h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$HOSPITALONSET)[2])
outofhospitalonset8h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$HOSPITALONSET)[2])
outofhospitalonset9h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$HOSPITALONSET)[2])
outofhospitalonset10h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$HOSPITALONSET)[2])
outofhospitalonset11h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$HOSPITALONSET)[2])
outofhospitalonset12h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$HOSPITALONSET)[2])
outofhospitalonset13h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$HOSPITALONSET)[2])
outofhospitalonset14h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$HOSPITALONSET)[2])
outofhospitalonset15h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$HOSPITALONSET)[2])
outofhospitalonset16h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$HOSPITALONSET)[2])
outofhospitalonset17h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$HOSPITALONSET)[2])
outofhospitalonset18h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$HOSPITALONSET)[2])
outofhospitalonset19h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$HOSPITALONSET)[2])
outofhospitalonset20h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$HOSPITALONSET)[2])
outofhospitalonset21h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$HOSPITALONSET)[2])
outofhospitalonset22h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$HOSPITALONSET)[2])
outofhospitalonset23h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$HOSPITALONSET)[2])
outofhospitalonset24h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$HOSPITALONSET)[2])


circadiandata$outofhospitalonset <- c(outofhospitalonset1h, outofhospitalonset2h, outofhospitalonset3h,
                                      outofhospitalonset4h, outofhospitalonset5h, outofhospitalonset6h,
                                      outofhospitalonset7h, outofhospitalonset8h, outofhospitalonset9h,
                                      outofhospitalonset10h, outofhospitalonset11h, outofhospitalonset12h,
                                      outofhospitalonset13h, outofhospitalonset14h, outofhospitalonset15h,
                                      outofhospitalonset16h, outofhospitalonset17h, outofhospitalonset18h,
                                      outofhospitalonset19h, outofhospitalonset20h, outofhospitalonset21h,
                                      outofhospitalonset22h, outofhospitalonset23h, outofhospitalonset24h)
circadiandata$outofhospitalonset <- ifelse(is.na(circadiandata$outofhospitalonset) == TRUE, 0, circadiandata$outofhospitalonset)
circadiandata$numberoutofhospitalonset <- circadiandata$numberSE * circadiandata$outofhospitalonset

## Cosinor analysis
SEfitcomparisonhospitalonset <- cosinor.lm(numberSE ~ time(hourofSE) + numberoutofhospitalonset + amp.acro(numberoutofhospitalonset), data = circadiandata, period = 24)
summary(SEfitcomparisonhospitalonset)

## Raw model coefficients:
##                              estimate standard.error lower.CI upper.CI
## (Intercept)                    6.9350         3.4894   0.0958  13.7741
## numberoutofhospitalonset       0.7835         0.3002   0.1950   1.3719
## rrr                           -1.6279         4.2202  -9.8993   6.6436
## sss                           -2.3001         4.4855 -11.0916   6.4914
## numberoutofhospitalonset:rrr   0.0168         0.3766  -0.7213   0.7549
## numberoutofhospitalonset:sss   0.2380         0.3821  -0.5109   0.9869
##                              p.value
## (Intercept)                   0.0469
## numberoutofhospitalonset      0.0091
## rrr                           0.6997
## sss                           0.6081
## numberoutofhospitalonset:rrr  0.9645
## numberoutofhospitalonset:sss  0.5334
## 
## ***********************
## 
## Transformed coefficients:
##                                    estimate standard.error lower.CI
## (Intercept)                          6.9350         3.4894   0.0958
## [numberoutofhospitalonset = 1]       0.7835         0.3002   0.1950
## amp                                  2.8179         4.1158  -5.2490
## [numberoutofhospitalonset = 1]:amp   2.6169         3.7532  -4.7393
## acr                                  0.9549         1.6259  -2.2317
## [numberoutofhospitalonset = 1]:acr   0.9076         1.6069  -2.2419
##                                    upper.CI p.value
## (Intercept)                         13.7741  0.0469
## [numberoutofhospitalonset = 1]       1.3719  0.0091
## amp                                 10.8848  0.4936
## [numberoutofhospitalonset = 1]:amp   9.9730  0.4857
## acr                                  4.1415  0.5570
## [numberoutofhospitalonset = 1]:acr   4.0570  0.5722

test_cosinor(SEfitcomparisonhospitalonset, "numberoutofhospitalonset", param = "amp")

## Global test: 
## Statistic: 
## [1] 0.23
## 
## 
##  P-value: 
## [1] 0.6307
## 
##  Individual tests: 
## Statistic: 
## [1] -0.48
## 
## 
##  P-value: 
## [1] 0.6307
## 
##  Estimate and confidence interval[1] "-0.2 (-1.02 to 0.62)"

test_cosinor(SEfitcomparisonhospitalonset, "numberoutofhospitalonset", param = "acr")

## Global test: 
## Statistic: 
## [1] 0.39
## 
## 
##  P-value: 
## [1] 0.5311
## 
##  Individual tests: 
## Statistic: 
## [1] -0.63
## 
## 
##  P-value: 
## [1] 0.5311
## 
##  Estimate and confidence interval[1] "-0.05 (-0.2 to 0.1)"

##Graph
hospitalonsetfigure <- ggplot.cosinor.lm(SEfitcomparisonhospitalonset, x_str = "numberoutofhospitalonset") +
  geom_bar(data = pSERGcircadian[pSERGcircadian$HOSPITALONSET == "no", ], aes(x = hourofSE, y = events), stat = 'identity', color = "#5698a3", fill = "#5698a3", alpha = 0.2) + 
  geom_bar(data = pSERGcircadian[pSERGcircadian$HOSPITALONSET == "yes", ], aes(x = hourofSE, y = events), stat = 'identity', color = "#a30234", fill = "#a30234", alpha = 0.2) +
    scale_x_discrete(limits = c(1,24)) +
  theme(panel.background=element_rect(fill = "white"), 
        axis.text=element_text(size = 18, color = "black", face = "bold"),
        axis.title=element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length=unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
        legend.position = "none") +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonhospitalonset$coefficients[1], yend = SEfitcomparisonhospitalonset$coefficients[1]), color = "#a30234") +
  geom_segment(aes(x = 15.7, xend = 15.7, y = SEfitcomparisonhospitalonset$coefficients[1], yend = (SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[3])), color = "#a30234", size=1.2) +
  geom_segment(aes(x = 3.7, xend = 3.7, y = SEfitcomparisonhospitalonset$coefficients[1], yend = (SEfitcomparisonhospitalonset$coefficients[1] - SEfitcomparisonhospitalonset$coefficients[3])), color = "#a30234", size=1.2) +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2], yend = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2]), color = "#5698a3") +
  geom_segment(aes(x = 15.4, xend = 15.4, y = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2], yend = (SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2] + SEfitcomparisonhospitalonset$coefficients[4])), color = "#5698a3", size=1.2) + geom_segment(aes(x = 3.4, xend = 3.4, y = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2], yend = (SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2] - SEfitcomparisonhospitalonset$coefficients[4])), color = "#5698a3", size=1.2)

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

hospitalonsetfigure

# Interactive graph
ggplotly(hospitalonsetfigure)

Time to first BZD

##Transform hours of SE to factor
pSERGcircadian$hourofSE<- as.factor(pSERGcircadian$hourofSE)

# Transform time to first BZD time to numeric
pSERGcircadian$BZDTIME.0 <- as.numeric(as.character(pSERGcircadian$BZDTIME.0))

## Warning: NAs introduced by coercion

#Summary data
summary(pSERGcircadian$BZDTIME.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00    5.00   16.50   70.06   45.00 1605.00      28

######################Transform hours of SE to numeric############
pSERGcircadian$hourofSE <- as.numeric(pSERGcircadian$hourofSE)

#Cosinor analysis
BZDfit <- cosinor.lm(BZDTIME.0 ~ time(hourofSE), data = pSERGcircadian[!is.na(pSERGcircadian$hourofSE), ], period = 24)
summary(BZDfit)

## Raw model coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  71.1615        10.2861  51.0012  91.3218  0.0000
## rrr          10.3922        14.6235 -18.2693  39.0536  0.4773
## sss           0.4276        14.3476 -27.6932  28.5483  0.9762
## 
## ***********************
## 
## Transformed coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  71.1615        10.2861  51.0012  91.3218  0.0000
## amp          10.4010        14.6555 -18.3233  39.1252  0.4779
## acr           0.0411         1.3763  -2.6564   2.7386  0.9762

##Draw the number of SE in a plot
################Transform hours of SE to factor######
pSERGcircadian$hourofSE <- as.factor(pSERGcircadian$hourofSE)
timeBZD <- ggplot.cosinor.lm(BZDfit) +
  geom_boxplot(data = pSERGcircadian[!is.na(pSERGcircadian$hourofSE), ], aes(x = hourofSE, y = BZDTIME.0), color = "#0076c0", fill = "#67771a", alpha = 0.5) + 
  scale_x_discrete(limits = c(1,24)) + 
  theme(panel.background=element_rect(fill = "white"), 
        axis.text=element_text(size = 16, color = "black", face = "bold"),
        axis.title=element_text(size = 18, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_discrete(breaks = c("4", "8", "12", "16", "20", "24")) + 
  scale_y_continuous(breaks = c(0, 50, 100, 150, 200, 250, 300, 350)) + 
  coord_cartesian(ylim = c(0, 200)) +
  labs(x = "Time of day (hour)", y= "Time to the first BZD (min)") +
  geom_segment(aes(x = 0,xend = 24, y= BZDfit$coefficients[1], yend = BZDfit$coefficients[1])) +
  geom_segment(aes(x = 0.5, xend = 0.5, y = BZDfit$coefficients[1], yend = BZDfit$coefficients[1] + BZDfit$coefficients[2]), size=1.2) +
  geom_segment(aes(x = 12.5, xend = 12.5, y = BZDfit$coefficients[1], yend = BZDfit$coefficients[1] - BZDfit$coefficients[2]), size=1.2)

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

timeBZD

## Warning: Removed 28 rows containing non-finite values (stat_boxplot).

# Interactive graph
ggplotly(timeBZD)

## Warning: Removed 28 rows containing non-finite values (stat_boxplot).

Time to first non-BZD AED

## Transform hours of SE to factor
pSERGcircadian$hourofSE <- as.factor(pSERGcircadian$hourofSE)


# Transform time to first AED time to numeric
pSERGcircadian$AEDTIME.0 <- as.numeric(as.character(pSERGcircadian$AEDTIME.0))

## Warning: NAs introduced by coercion

# Summary data
summary(pSERGcircadian$AEDTIME.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0    35.0    66.0   165.4   150.0  4320.0      27

## Transform hours of SE to numeric 
pSERGcircadian$hourofSE <- as.numeric(pSERGcircadian$hourofSE)

# Cosinor analysis
ASMfit <- cosinor.lm(AEDTIME.0~time(hourofSE), data=pSERGcircadian, period=24)
summary(ASMfit)

## Raw model coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept) 171.9320        19.5132 133.6869 210.1771  0.0000
## rrr          50.7544        28.0200  -4.1639 105.6727  0.0701
## sss          44.6559        26.9189  -8.1042  97.4159  0.0971
## 
## ***********************
## 
## Transformed coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept) 171.9320        19.5132 133.6869 210.1771  0.0000
## amp          67.6029        28.5555  11.6352 123.5706  0.0179
## acr           0.7216         0.3898  -0.0424   1.4855  0.0641

################Transform hours of SE to factor######
pSERGcircadian$hourofSE <- as.factor(pSERGcircadian$hourofSE)

#Draw the data and the fitting curve
timeASM <- ggplot.cosinor.lm(ASMfit) +
  geom_boxplot(data = pSERGcircadian[!is.na(pSERGcircadian$hourofSE), ], aes(hourofSE, AEDTIME.0), color = "#0076c0", fill = "#67771a", alpha = 0.5) +
  scale_x_discrete(limits = c(1,24)) + 
  scale_y_discrete(limits = c(0,800)) + 
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_text(size = 16, color = "black", face = "bold"),
        axis.title = element_text(size = 18, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_discrete(breaks = c("4", "8", "12", "16", "20", "24")) +
  scale_y_continuous(breaks = c(0, 100, 200, 300, 400, 500, 600)) + 
  coord_cartesian(ylim = c(0, 600)) +
  labs(x = "Time of day (hour)", y = "Time to the first non-BZD AED (min)") +
  geom_segment(aes(x = 0, xend = 24, y = ASMfit$coefficients[1], yend = ASMfit$coefficients[1])) +
  geom_segment(aes(x = 2.75, xend = 2.75, y = ASMfit$coefficients[1], yend = ASMfit$coefficients[1] + ASMfit$coefficients[2]), size=1.2) +
  geom_segment(aes(x = 14.75, xend = 14.75, y = ASMfit$coefficients[1], yend = ASMfit$coefficients[1] - ASMfit$coefficients[2]), size=1.2)

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

timeASM

## Warning: Removed 27 rows containing non-finite values (stat_boxplot).

# Interactive graph
ggplotly(timeASM)

## Warning: Removed 27 rows containing non-finite values (stat_boxplot).

Alternative hour definition

#################Create variable sensitivity analysis hour of SE######################
## Define the newly defined hour
pSERGcircadian$minutes <- pSERGcircadian$TIMESEIZURE_MINUTES
pSERGcircadian$minutesashours <- pSERGcircadian$minutes / 60
pSERGcircadian$minutesashours <- ifelse(is.na(pSERGcircadian$minutesashours) == TRUE | pSERGcircadian$minutesashours > 1, 0, pSERGcircadian$minutesashours)

## Create variable new hour
pSERGcircadian$newhour <- as.numeric(pSERGcircadian$hourofSE) + pSERGcircadian$minutesashours

## Rename hours
pSERGcircadian$newh[pSERGcircadian$newhour >= 24.5 | pSERGcircadian$newhour < 1.5] <- 1
pSERGcircadian$newh[pSERGcircadian$newhour >= 1.5 & pSERGcircadian$newhour < 2.5] <- 2
pSERGcircadian$newh[pSERGcircadian$newhour >= 2.5 & pSERGcircadian$newhour < 3.5] <- 3
pSERGcircadian$newh[pSERGcircadian$newhour >= 3.5 & pSERGcircadian$newhour < 4.5] <- 4
pSERGcircadian$newh[pSERGcircadian$newhour >= 4.5 & pSERGcircadian$newhour < 5.5] <- 5
pSERGcircadian$newh[pSERGcircadian$newhour >= 5.5 & pSERGcircadian$newhour < 6.5] <- 6
pSERGcircadian$newh[pSERGcircadian$newhour >= 6.5 & pSERGcircadian$newhour < 7.5] <- 7
pSERGcircadian$newh[pSERGcircadian$newhour >= 7.5 & pSERGcircadian$newhour < 8.5] <- 8
pSERGcircadian$newh[pSERGcircadian$newhour >= 8.5 & pSERGcircadian$newhour < 9.5] <- 9
pSERGcircadian$newh[pSERGcircadian$newhour >= 9.5 & pSERGcircadian$newhour < 10.5] <- 10
pSERGcircadian$newh[pSERGcircadian$newhour >= 10.5 & pSERGcircadian$newhour < 11.5] <- 11
pSERGcircadian$newh[pSERGcircadian$newhour >= 11.5 & pSERGcircadian$newhour < 12.5] <- 12
pSERGcircadian$newh[pSERGcircadian$newhour >= 12.5 & pSERGcircadian$newhour < 13.5] <- 13
pSERGcircadian$newh[pSERGcircadian$newhour >= 13.5 & pSERGcircadian$newhour < 14.5] <- 14
pSERGcircadian$newh[pSERGcircadian$newhour >= 14.5 & pSERGcircadian$newhour < 15.5] <- 15
pSERGcircadian$newh[pSERGcircadian$newhour >= 15.5 & pSERGcircadian$newhour < 16.5] <- 16
pSERGcircadian$newh[pSERGcircadian$newhour >= 16.5 & pSERGcircadian$newhour < 17.5] <- 17
pSERGcircadian$newh[pSERGcircadian$newhour >= 17.5 & pSERGcircadian$newhour < 18.5] <- 18
pSERGcircadian$newh[pSERGcircadian$newhour >= 18.5 & pSERGcircadian$newhour < 19.5] <- 19
pSERGcircadian$newh[pSERGcircadian$newhour >= 19.5 & pSERGcircadian$newhour < 20.5] <- 20
pSERGcircadian$newh[pSERGcircadian$newhour >= 20.5 & pSERGcircadian$newhour < 21.5] <- 21
pSERGcircadian$newh[pSERGcircadian$newhour >= 21.5 & pSERGcircadian$newhour < 22.5] <- 22
pSERGcircadian$newh[pSERGcircadian$newhour >= 22.5 & pSERGcircadian$newhour < 23.5] <- 23
pSERGcircadian$newh[pSERGcircadian$newhour >= 23.5 & pSERGcircadian$newhour < 24.5] <- 24


# Test (Kolmogorov-Smirnov) if the distribution of SE during the day is consistent with a uniform distribution
ks.test(pSERGcircadian$newh, "punif", 0, 24)

## Warning in ks.test(pSERGcircadian$newh, "punif", 0, 24): ties should not be
## present for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  pSERGcircadian$newh
## D = 0.082428, p-value = 0.01347
## alternative hypothesis: two-sided

#####################################################################################
# Create another database with the circadian data
newoccurrences <- count(pSERGcircadian, "newh")
# Make it a dataframe
newcircadiandata <- data.frame(newoccurrences)
newcircadiandata$numberSE <- newcircadiandata$freq
newcircadiandata <- newcircadiandata[ , c("newh", "numberSE")]
#####################################################################################



# Cosinor analysis
newSEfit <- cosinor.lm(numberSE ~ time(newh), data = newcircadiandata, period = 24)
summary(newSEfit)

## Raw model coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  15.3333         0.7375  13.8878  16.7788  0.0000
## rrr          -3.4051         1.0430  -5.4493  -1.3608  0.0011
## sss          -0.0730         1.0430  -2.1172   1.9712  0.9442
## 
## ***********************
## 
## Transformed coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  15.3333         0.7375  13.8878  16.7788  0.0000
## amp           3.4058         1.0430   1.3616   5.4501  0.0011
## acr           0.0214         0.3062  -0.5788   0.6216  0.9442

################Every hour new definition ###############
## Draw the number of SE episodes in a plot with stacked units
newdistributionhour <- ggplot(pSERGcircadian, aes(x = round(newh), y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") + 
  scale_x_discrete(limits = c(1,24)) +
  coord_cartesian(xlim = c(1, 24)) +
  theme(panel.background = element_rect(fill = "white"), 
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15, 20, 25)) +
  labs(x= "Time of day (hour)", y= "Number of rSE episodes")

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

newdistributionhour

# Interactive graph
ggplotly(newdistributionhour)

# Draw the fitting curve 1h
## Draw the number of SE in a plot with the redefined hour
newcircadianfigure <- ggplot.cosinor.lm(SEfit) + geom_bar(data = pSERGcircadian, aes(x = newh, y = events), stat = 'identity', color = "#0076c0", fill = "#67771a", alpha = 0.8) + 
  scale_x_discrete(limits = c(1,24)) +
  theme(panel.background = element_rect(fill = "white"), 
        axis.text=element_text(size = 18, color = "black", face = "bold"),
        axis.title=element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15, 20, 25)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
  geom_segment(aes(x = 0, xend = 24, y = SEfit$coefficients[1], yend = SEfit$coefficients[1])) +
  geom_segment(aes(x = 11.7, xend = 11.7, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] + SEfit$coefficients[2])), size=1.2) +
  geom_segment(aes(x = 23.7, xend = 23.7, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] - SEfit$coefficients[2])), size=1.2)

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

newcircadianfigure

# Interactive graph
ggplotly(newcircadianfigure)